#include "bc_orthogonalFlowVelocity_waterdepth.hpp"
#include "meshTypes.hpp"
#include "field.hpp"
#include "boundaryPatch.hpp"
#include "sweFunctions.hpp"
#include "simulationTime.hpp"
#include "exceptions.hpp"
#include "geometry.hpp"

namespace hms {

#define CHEAD(RETURN_TYPE) \
template<typename MeshType> \
RETURN_TYPE OrthogonalFlowVelocity_waterdepth<MeshType>::

CHEAD() OrthogonalFlowVelocity_waterdepth(
	Field<MeshType>& f,
	const BoundaryPatch& p,
	const Field<MeshType>& vFieldArg,
	scalar vBoundaryArg
) :
	BaseType(Type::orthogonalFlowVelocity_waterdepth, f, p)
{
	vField(vFieldArg);
	vBoundary(vBoundaryArg);
}

CHEAD(void) update() const {
	Field<MeshType>& d { field() };
	const Field<MeshType>& v { vField() };
	const BoundaryPatch& patch { this->patch() };
	/* Simons (2019) states that for supercritical flow, dGhost has to be set
	 * to dInner (effectively zeroGradient). However, in the accompanying 
	 * software hms, this is not implemented. The implementation of hms is
	 * mimicked here. */
	
	/* in Simons (2019), vBoundary has the opposite sign (positive for
	 * outflow, negative for inflow). Therefore, the sign of m_vBoundary is
	 * inverted here in comparison to Simons (2019), Eq. 3.133. */
	auto setValues = [&]( const auto& vOrthoDomain ){
		d.ghostValues(patch) = (
				vOrthoDomain +
				vBoundary() +
				2 * sqrt(phys::g * d.domainValues(patch))
			).square() / 
			(4*phys::g);
	};

	/* rectangular meshes
	 * (the non-rectangular approach works for rectangular meshes as well,
	 * but is a bit slower) */
	if constexpr (
		std::is_same_v<MeshType, UniMesh> ||
		std::is_same_v<MeshType, RectMesh>
	){
		using BL = BoundaryLocation;
		BL location { whereIs( patch, v.mesh().size() ) };

		int
			orthoRow { location == BL::bottom || location == BL::top ?  1 : 0 },
			sign     { location == BL::bottom || location == BL::left ? -1 : 1 };
		
		setValues( sign * v.domainValues(patch).row(orthoRow) );

	/* non-rectangular meshes (Mesh, StructMesh) */
	} else {
		setValues(
			( // projection of v onto n is edge orthogonal component vOrthoDomain
				v.domainValues(patch) *
				v.mesh().edges().outwardNormals(patch).array()
			).colwise().sum()
		);
	}
}

CHEAD(auto) vBoundary() const -> scalar {
	return m_vBoundary;
}

CHEAD(void) vBoundary(scalar vBoundaryArg){
	m_vBoundary = vBoundaryArg;
}

#undef CHEAD

#define INSTANTIATE_FOR_ALL_MESHTYPES(MeshType) \
template class OrthogonalFlowVelocity_waterdepth<MeshType>;
#include "meshInstantiator.hpp"

} // end namespace hms